#############################################################################
# Figure2_FigureB1_FigureB2.R
# Aim: to replicate Figure 2, Figure B1, and Figure B2 in Atsusaka and Stevenson (2021)
# Run time: about 47 minutes (ALARM: LONG RUN TIME!)
#############################################################################

library(tidyverse)
rm(list=ls())

start <- Sys.time()

Nsim <- 2000         # THIS MUST BE 8000 IN THE FINAL VERSION
naive.cover <- NA
bc.cover <- NA
enzmann.cover1 <- NA
enzmann.cover2 <- NA
enzmann.cover3 <- NA
schnapp.cover1 <- NA
schnapp.cover2 <- NA
schnapp.cover3 <- NA
naive.pred <- NA
bc.pred <- NA
enzmann.pred1 <- NA
enzmann.pred2 <- NA
enzmann.pred3 <- NA
schnapp.pred1 <- NA
schnapp.pred2 <- NA
schnapp.pred3 <- NA

pi.true <- NA
Sim.N <- NA
Sim.pi <- NA
Sim.p <- NA
Sim.p2 <- NA
Sim.gamma <- NA 
REbc <- NA
REenzmann <- NA

# STORAGES OF RESULTS (k)
RElist.bc <- list()
RElist.enzmann <- list()
sim.bias.naive <- numeric()
sim.bias.bc <- numeric()
sim.bias.enzmann1 <- numeric()
sim.bias.enzmann2 <- numeric()
sim.bias.enzmann3 <- numeric()
sim.bias.schnapp1 <- numeric()
sim.bias.schnapp2 <- numeric()
sim.bias.schnapp3 <- numeric()
sim.rmse.naive <- numeric()
sim.rmse.bc <- numeric()
sim.rmse.enzmann1 <- numeric()
sim.rmse.enzmann2 <- numeric()
sim.rmse.enzmann3 <- numeric()
sim.rmse.schnapp1 <- numeric()
sim.rmse.schnapp2 <- numeric()
sim.rmse.schnapp3 <- numeric()
sim.cover.naive <- numeric()
sim.cover.bc <- numeric()
sim.cover.enzmann1 <- numeric()
sim.cover.enzmann2 <- numeric()
sim.cover.enzmann3 <- numeric()
sim.cover.schnapp1 <- numeric()
sim.cover.schnapp2 <- numeric()
sim.cover.schnapp3 <- numeric()

pb <- txtProgressBar(min=1, max=Nsim, initial=0, style=3) # PROGRESS BAR

sample.size = c(200, 500, 1000, 2000, 5000)
set.seed(12345)
for(k in 1:5){
N=sample.size[k] # CHOOSE THE SAMPLE SIZE

for(i in 1:Nsim){

setTxtProgressBar(pb,i)

  # (1) DATA GENERATING PROCESS
  pi = runif(1, min=0.01, max=0.45)           # Plausible true prevelance
  p  = runif(1, min=0.088, max=0.33)          # Plausible randomization probability
  p2 = runif(1, min=0.088, max=0.33)          # Plausible randomization probability for anchor question
  gamma = runif(1, min=0.5, max=1)            # Plausible attentive rate (only consider over 0.5)
  Attentive = rbinom(n=N, size=1, prob=gamma) # Attentive R or not
  
# SENSITIVE QUESTION OF INTERST
  StatementA = rbinom(n=N, size=1, prob=pi)           # Sesitive item
  StatementB = rbinom(n=N, size=1, prob=p)            # Randomization item
  True.res = ifelse(StatementA != StatementB, 0, 1)   # Survey answer 
  
  Obs.res = rep(NA, N)
  Obs.res[Attentive==1] = True.res[Attentive==1]      # Observed answer for attentive Rs
  Obs.res[Attentive==0] = rbinom(n=length(N[Attentive==0]),
                                 size=1, prob=0.5)           # Observed answer, otherwise
  
# NON-SENSITIVE AUXIALLY QUESTION
  StatementC = 0                                      # Attention "c"heck item (True prevalence is 0 for everyone)
  StatementD = rbinom(n=N, size=1, prob=p2)           # Anchor question
  True.res.prime = ifelse(StatementC != StatementD, 0, 1)   # Survey answer in Anchor Question
  
  Obs.res.prime = rep(NA, N)
  Obs.res.prime[Attentive==1] = True.res.prime[Attentive==1] # Observed answer for attentive Rs
  Obs.res.prime[Attentive==0] = rbinom(n=length(N[Attentive==0]),
                                       size=1, prob=0.5)            # Observed answer, otherwise
  
  dat <- cbind(Obs.res, True.res, Attentive, Obs.res.prime, True.res.prime)
  dat <- as_tibble(dat)
  head(dat) # See the D
  
  
# SCHNAPP-DIRECT QUESTION ("Did you give random answers in the last question?")
# Assumption S1. One Side Lying (Attentive Respondents Do Not Lie)
#                                Only Inattentive Respondents Lie at Some Rate

# Lying Parameter 
  lying <- c(0, 0.25, 0.5) # 10%, 25%, 50% are lying among inattentive respondents
   
  SDQ1 <- NA
  SDQ1[Attentive==1] <- 0 # 0 IF NOT GIVING RANDOM ANSWERS
  SDQ1[Attentive==0] <- 1 # 1 IF GIVING RANDOM ANSWERS
  
  SDQ2 <- NA
  SDQ2[Attentive==1] <- 0 # 0 IF NOT GIVING RANDOM ANSWERS
  SDQ2[Attentive==0] <- rbinom(n=length(N[Attentive==0]), size=1, prob=1-lying[2])
  
  SDQ3 <- NA
  SDQ3[Attentive==1] <- 0
  SDQ3[Attentive==0] <- rbinom(n=length(N[Attentive==0]), size=1, prob=1-lying[3])  

  
# (2) NAIVE CROSSWISE MODEL
  lambda.hat = mean(Obs.res) # Observed proportion of YESYES or NONO
  pi.hat.naive = (lambda.hat+p-1)/(2*p-1)
  
  pi.hat.naive.var = (pi.hat.naive*(1-pi.hat.naive))/(N-1) +
    (p*(1-p))/((N-1)*((2*p-1)^2))
  pi.hat.naive.sd = sqrt(pi.hat.naive.var)     
  
  naive.low  = pi.hat.naive-1.96*pi.hat.naive.sd; naive.low  = ifelse(naive.low > 0,  naive.low, 0)
  naive.high = pi.hat.naive+1.96*pi.hat.naive.sd; naive.high = ifelse(naive.high < 1, naive.high, 1)
  
  
# (3) BIAS CORRECTED CROSSWISE MODEL
  gamma.hat = (mean(Obs.res.prime)-0.5)/(0.5-p2) # Estimated level of inattentiveness
  Bias.hat = (1/2)*((lambda.hat-0.5)/(p-0.5)) - (1/(2*gamma.hat))*((lambda.hat-0.5)/(p-0.5))
  Bias.hat
  
  pi.hat.bc = pi.hat.naive - Bias.hat # Bias Correction
  pi.hat.bc <- ifelse(pi.hat.bc < 0, 0, ifelse(pi.hat.bc > 1, 1, pi.hat.bc))
  
# BOOTSTRAPPING (1000 TIMES) FOR INFERENCE
  bs <- NA
  for(b in 1:1000){
    index <- sample(1:nrow(dat), size=dim(dat), replace=TRUE)  
    bs.dat <- dat[index,]                                      
    
    bs.lambda.hat = mean(bs.dat$Obs.res)                    # Observed proportion of TRUE-TRUE or FALSE-FALSE
    bs.pi.hat.naive = (bs.lambda.hat+p-1)/(2*p-1)
    bs.gamma.hat = (mean(bs.dat$Obs.res.prime)-0.5)/(0.5-p2) # Estimated level of inattentiveness
    bs.bias.hat = (1/2)*((bs.lambda.hat-0.5)/(p-0.5)) - (1/(2*bs.gamma.hat))*((bs.lambda.hat-0.5)/(p-0.5)) 
    
    bs[b] = bs.pi.hat.naive - bs.bias.hat # Bias Correction within Bootstrapping
    
    bs[b] <- ifelse(bs[b] < 0, 0, bs[b])
    bs[b] <- ifelse(bs[b] > 1, 1, bs[b])
  }
  
  pi.hat.bc.var = var(bs)
  pi.hat.bc.sd = sqrt(pi.hat.bc.var)
  
  bc.low  = quantile(bs, prob=0.025, na.rm=T); bc.low  = ifelse(bc.low > 0,  bc.low, 0)
  bc.high = quantile(bs, prob=0.975, na.rm=T); bc.high = ifelse(bc.high < 1, bc.high, 1)
  
  
# (4) ENZMANN-SCHNAPP CORRECTION WITH THREE LYING PARAMETERS
# 0%  Liars
r = mean(SDQ1) # ESTIMATED PROPORTION OF RANDOM ANSWERS VIA DIRECT QUESTIONING
pi.hat.enzmann1 = (pi.hat.naive - 0.5*r)/(1-r)

A = pi.hat.naive*(1-pi.hat.naive)/(N-1)
B = p*(1-p)/ ((N-1)*((2*p - 1)^2))
C = (1-gamma.hat)^2

pi.hat.enzmann.var = (A+B)/C
pi.hat.enzmann.sd = sqrt(pi.hat.enzmann.var)

enzmann.low1  = pi.hat.enzmann1-1.96*pi.hat.enzmann.sd
enzmann.high1 = pi.hat.enzmann1+1.96*pi.hat.enzmann.sd 


# 25%  Liars
r = mean(SDQ2) # ESTIMATED PROPORTION OF RANDOM ANSWERS VIA DIRECT QUESTIONING
pi.hat.enzmann2 = (pi.hat.naive - 0.5*r)/(1-r)

A = pi.hat.naive*(1-pi.hat.naive)/(N-1)
B = p*(1-p)/ ((N-1)*((2*p - 1)^2))
C = (1-gamma.hat)^2

pi.hat.enzmann.var = (A+B)/C
pi.hat.enzmann.sd = sqrt(pi.hat.enzmann.var)

enzmann.low2  = pi.hat.enzmann2-1.96*pi.hat.enzmann.sd
enzmann.high2 = pi.hat.enzmann2+1.96*pi.hat.enzmann.sd 


# 50%  Liars
r = mean(SDQ3) # ESTIMATED PROPORTION OF RANDOM ANSWERS VIA DIRECT QUESTIONING
pi.hat.enzmann3 = (pi.hat.naive - 0.5*r)/(1-r)

A = pi.hat.naive*(1-pi.hat.naive)/(N-1)
B = p*(1-p)/ ((N-1)*((2*p - 1)^2))
C = (1-r)^2

pi.hat.enzmann.var = (A+B)/C
pi.hat.enzmann.sd = sqrt(pi.hat.enzmann.var)

enzmann.low3  = pi.hat.enzmann3-1.96*pi.hat.enzmann.sd
enzmann.high3 = pi.hat.enzmann3+1.96*pi.hat.enzmann.sd



# (5) SCHNAPP CORRECTION (CMR-I)
# 0%, 50%, 90% are Lying
Adj.res1 <- ifelse(SDQ1==0, Obs.res, 0.5) # Otherwise, change the answer to 0.5
Adj.res2 <- ifelse(SDQ2==0, Obs.res, 0.5) # Otherwise, change the answer to 0.5
Adj.res3 <- ifelse(SDQ3==0, Obs.res, 0.5) # Otherwise, change the answer to 0.5

Adj.lambda1 <- mean(Adj.res1)
Adj.lambda2 <- mean(Adj.res2)
Adj.lambda3 <- mean(Adj.res3)

Adj.naive1 = (Adj.lambda1 +p-1)/(2*p-1)
Adj.naive2 = (Adj.lambda2 +p-1)/(2*p-1)
Adj.naive3 = (Adj.lambda3 +p-1)/(2*p-1)

pi.hat.schnapp.var1 = (Adj.naive1*(1-Adj.naive1))/(N-1) +
                     (p*(1-p))/((N-1)*((2*p-1)^2))
pi.hat.schnapp.sd1 = sqrt(pi.hat.schnapp.var1)  
pi.hat.schnapp.var2 = (Adj.naive2*(1-Adj.naive2))/(N-1) +
                     (p*(1-p))/((N-1)*((2*p-1)^2))
pi.hat.schnapp.sd2 = sqrt(pi.hat.schnapp.var2)  
pi.hat.schnapp.var3 = (Adj.naive3*(1-Adj.naive3))/(N-1) +
                     (p*(1-p))/((N-1)*((2*p-1)^2))
pi.hat.schnapp.sd3 = sqrt(pi.hat.schnapp.var3)  

schnapp.low1  = Adj.naive1-1.96*pi.hat.schnapp.sd1
schnapp.high1 = Adj.naive1+1.96*pi.hat.schnapp.sd1
schnapp.low2  = Adj.naive2-1.96*pi.hat.schnapp.sd2
schnapp.high2 = Adj.naive2+1.96*pi.hat.schnapp.sd2
schnapp.low3  = Adj.naive3-1.96*pi.hat.schnapp.sd3
schnapp.high3 = Adj.naive3+1.96*pi.hat.schnapp.sd3



# STORE OUTPUTS
  naive.pred[i] <- pi.hat.naive
  bc.pred[i] <- pi.hat.bc
  enzmann.pred1[i] <- pi.hat.enzmann1
  enzmann.pred2[i] <- pi.hat.enzmann2
  enzmann.pred3[i] <- pi.hat.enzmann3
  schnapp.pred1[i] <- Adj.naive1
  schnapp.pred2[i] <- Adj.naive2
  schnapp.pred3[i] <- Adj.naive3
  pi.true[i] <- pi
  naive.cover[i] <- (naive.low <= pi & pi <= naive.high) # 95CI contains pi
  bc.cover[i] <-    (bc.low <= pi & pi <= bc.high)       # 95CI contains pi
  enzmann.cover1[i]<-(enzmann.low1<=pi & pi<=enzmann.high1)  
  enzmann.cover2[i]<-(enzmann.low2<=pi & pi<=enzmann.high2)  
  enzmann.cover3[i]<-(enzmann.low3<=pi & pi<=enzmann.high3)  
  schnapp.cover1[i]<-(schnapp.low1<=pi & pi<=schnapp.high1)  
  schnapp.cover2[i]<-(schnapp.low2<=pi & pi<=schnapp.high2)  
  schnapp.cover3[i]<-(schnapp.low3<=pi & pi<=schnapp.high3)  
  
  Sim.N[i] <- N
  Sim.pi[i] <- pi
  Sim.p[i] <- p 
  Sim.p2[i] <- p2
  Sim.gamma[i] <- gamma 
  
  REbc[i] <- (bc.high-bc.low)/(naive.high-naive.low)                # LENGTH OF 95CI
  REenzmann[i] <- (enzmann.high1-enzmann.low1)/(naive.high-naive.low) # LENGTH OF 95CI
}

# Bias: sum_{i=1}^{n}(\hat\pi - \pi)/n
bias.naive.cm      <- sum(naive.pred - pi.true)/Nsim 
bias.bias.correct  <- sum(bc.pred - pi.true)/Nsim
bias.enzmann1       <- sum(enzmann.pred1 - pi.true)/Nsim
bias.enzmann2       <- sum(enzmann.pred2 - pi.true)/Nsim
bias.enzmann3       <- sum(enzmann.pred3 - pi.true)/Nsim
bias.schnapp1       <- sum(schnapp.pred1 - pi.true)/Nsim
bias.schnapp2       <- sum(schnapp.pred2 - pi.true)/Nsim
bias.schnapp3       <- sum(schnapp.pred3 - pi.true)/Nsim

# RMSE: \sqrt( (sum_{i=1}^{n}(\hat\pi - \pi)^2)/n )
rmse.naive.cm      <- sum( (naive.pred - pi.true)^2 )/Nsim
rmse.bias.correct  <- sum( (bc.pred - pi.true)^2 )/Nsim
rmse.enzmann1        <- sum( (enzmann.pred1 - pi.true)^2 )/Nsim
rmse.enzmann2        <- sum( (enzmann.pred1 - pi.true)^2 )/Nsim
rmse.enzmann3        <- sum( (enzmann.pred1 - pi.true)^2 )/Nsim
rmse.schnapp1       <- sum( (schnapp.pred1 - pi.true)^2 )/Nsim
rmse.schnapp2       <- sum( (schnapp.pred2 - pi.true)^2 )/Nsim
rmse.schnapp3       <- sum( (schnapp.pred3 - pi.true)^2 )/Nsim


# What Percentage of Times Did 95CIs Capture the Ground Truth?
coverage.naive.cm     <- mean(naive.cover, na.rm =T)*100
coverage.bias.correct <- mean(bc.cover, na.rm=T)*100
coverage.enzmann1      <- mean(enzmann.cover1, na.rm=T)*100
coverage.enzmann2      <- mean(enzmann.cover1, na.rm=T)*100
coverage.enzmann3      <- mean(enzmann.cover1, na.rm=T)*100
coverage.schnapp1     <- mean(schnapp.cover1, na.rm=T)*100
coverage.schnapp2     <- mean(schnapp.cover2, na.rm=T)*100
coverage.schnapp3     <- mean(schnapp.cover3, na.rm=T)*100


# STORE RESULTS
# Relative Lenght of CI
RElist.bc[[k]] <- REbc
RElist.enzmann[[k]] <- REenzmann

sim.bias.naive <- rbind(sim.bias.naive, bias.naive.cm) 
sim.bias.bc <- rbind(sim.bias.bc, bias.bias.correct) 
sim.bias.enzmann1 <- rbind(sim.bias.enzmann1, bias.enzmann1)
sim.bias.enzmann2 <- rbind(sim.bias.enzmann2, bias.enzmann2) 
sim.bias.enzmann3 <- rbind(sim.bias.enzmann3, bias.enzmann3) 
sim.bias.schnapp1 <- rbind(sim.bias.schnapp1, bias.schnapp1)
sim.bias.schnapp2 <- rbind(sim.bias.schnapp2, bias.schnapp2) 
sim.bias.schnapp3 <- rbind(sim.bias.schnapp3, bias.schnapp3) 

sim.rmse.naive <- rbind(sim.rmse.naive, rmse.naive.cm) 
sim.rmse.bc <- rbind(sim.rmse.bc, rmse.bias.correct) 
sim.rmse.enzmann1 <- rbind(sim.rmse.enzmann1, rmse.enzmann1 ) 
sim.rmse.enzmann2 <- rbind(sim.rmse.enzmann2, rmse.enzmann2 ) 
sim.rmse.enzmann3 <- rbind(sim.rmse.enzmann3, rmse.enzmann3 ) 

sim.rmse.schnapp1 <- rbind(sim.rmse.schnapp1, rmse.schnapp1 ) 
sim.rmse.schnapp2 <- rbind(sim.rmse.schnapp2, rmse.schnapp2 ) 
sim.rmse.schnapp3 <- rbind(sim.rmse.schnapp3, rmse.schnapp3 ) 

sim.cover.naive <- rbind(sim.cover.naive, coverage.naive.cm) 
sim.cover.bc <- rbind(sim.cover.bc, coverage.bias.correct) 
sim.cover.enzmann1 <- rbind(sim.cover.enzmann1, coverage.enzmann1)
sim.cover.enzmann2 <- rbind(sim.cover.enzmann2, coverage.enzmann2)
sim.cover.enzmann3 <- rbind(sim.cover.enzmann3, coverage.enzmann3)
sim.cover.schnapp1 <- rbind(sim.cover.schnapp1, coverage.schnapp1)
sim.cover.schnapp2 <- rbind(sim.cover.schnapp2, coverage.schnapp2)
sim.cover.schnapp3 <- rbind(sim.cover.schnapp3, coverage.schnapp3)

} # END OF THE OUTER LOOP (k)


#############################################################################
# VISUALIZATION OF THE RESULTS
#############################################################################

# TUTN INTO VECTORS
sim.bias.naive <- sim.bias.naive %>% as_tibble() %>% pull() 
sim.bias.bc <- sim.bias.bc %>% as_tibble() %>% pull()
sim.bias.enzmann1 <- sim.bias.enzmann1 %>% as_tibble()%>% pull()
sim.bias.enzmann2 <- sim.bias.enzmann2 %>% as_tibble()%>% pull()
sim.bias.enzmann3 <- sim.bias.enzmann3 %>% as_tibble()%>% pull()
sim.bias.schnapp1 <- sim.bias.schnapp1 %>% as_tibble()%>% pull()
sim.bias.schnapp2 <- sim.bias.schnapp2 %>% as_tibble()%>% pull()
sim.bias.schnapp3 <- sim.bias.schnapp3 %>% as_tibble()%>% pull()
sim.rmse.naive <- sim.rmse.naive %>% as_tibble()%>% pull()
sim.rmse.bc <- sim.rmse.bc %>% as_tibble()%>% pull()
sim.rmse.enzmann1 <- sim.rmse.enzmann1 %>% as_tibble()%>% pull()
sim.rmse.enzmann2 <- sim.rmse.enzmann2 %>% as_tibble()%>% pull()
sim.rmse.enzmann3 <- sim.rmse.enzmann3 %>% as_tibble()%>% pull()
sim.rmse.schnapp1 <- sim.rmse.schnapp1 %>% as_tibble()%>% pull()
sim.rmse.schnapp2 <- sim.rmse.schnapp2 %>% as_tibble()%>% pull()
sim.rmse.schnapp3 <- sim.rmse.schnapp3 %>% as_tibble()%>% pull()
sim.cover.naive <- sim.cover.naive %>% as_tibble() %>% pull()
sim.cover.bc <- sim.cover.bc %>% as_tibble() %>% pull()
sim.cover.enzmann1 <- sim.cover.enzmann1 %>% as_tibble() %>% pull()
sim.cover.enzmann2 <- sim.cover.enzmann2 %>% as_tibble() %>% pull()
sim.cover.enzmann3 <- sim.cover.enzmann3 %>% as_tibble() %>% pull()
sim.cover.schnapp1 <- sim.cover.schnapp1 %>% as_tibble() %>% pull()
sim.cover.schnapp2 <- sim.cover.schnapp2 %>% as_tibble() %>% pull()
sim.cover.schnapp3 <- sim.cover.schnapp3 %>% as_tibble() %>% pull()


# VISUALIZATION
######################################################
# FIGURE 2
######################################################
pdf("Figure2.pdf", width=7.8, height=3)
par(mfrow=c(1,3), mar=c(2,3.5,4,0.5), oma=c(1.5,1,0,0.5)+0.02, mgp=c(3,0.6,0))
  
  plot(sim.bias.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, xlim=c(200,5000), ylim=c(-0.02, 0.08))
  abline(h=0, lty=2, col="firebrick4", lwd=2)
  lines(sim.bias.naive ~ sample.size, type="b", cex=1.8)
  lines(sim.bias.bc ~ sample.size, type="b", pch=19, cex=1.8)
  text(x=600, y=0.06, labels="Naive")
  text(x=1400, y=-0.015, labels="BC (Bias-Corrected)")
  title("Bias", line=1, cex.main=1.5)
  
  plot(sim.rmse.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, xlim=c(200,5000), ylim=c(0,0.012))
  abline(h=0, lty=2, col="firebrick4", lwd=2)
  lines(sim.rmse.naive ~ sample.size, type="b", cex=1.8)
  lines(sim.rmse.bc ~ sample.size, type="b", pch=19, cex=1.8)
  text(x=800, y=0.007, labels="Naive")
  text(x=500, y=0.0005, labels="BC")
  title("Root Mean Square Error", line=1, cex.main=1.5)
  
  
  plot(sim.cover.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, 
       xlim=c(200,5000), ylim=c(0,100))
  abline(h=95, lty=2, col="firebrick4", lwd=2)
  lines(sim.cover.naive ~ sample.size, type="b", cex=1.8)
  lines(sim.cover.bc ~ sample.size, type="b", pch=19, cex=1.8)
  text(x=600, y=35, labels="Naive")
  text(x=500, y=85, labels="BC")
  title("Coverage of 95%\nConfidence Interval", line=1, cex.main=1.5)
  
mtext("Sample Size", side=1, line=0, outer=T, cex=1)  
dev.off()
 

######################################################
# FIGURE B1
######################################################
pdf("FigureB1.pdf", width=7.8, height=3.2)
par(mfrow=c(1,3), mar=c(2,3.5,4,0.5), oma=c(1.5,1,0,0.5)+0.02, mgp=c(3,0.6,0))
  

col_vec <- c("black", "firebrick4",
             "seagreen2", "seagreen3", "seagreen4",
             "steelblue2","steelblue3","steelblue4")
pch_vec <- c(1,19,2,0,5,3,4,8)

  plot(sim.bias.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, xlim=c(200,5000), ylim=c(-0.12, 0.08))
  abline(h=0, lty=2, col="gray50", lwd=2)
  lines(sim.bias.naive ~ sample.size, type="b", cex=1.8, col=col_vec[1])
  lines(sim.bias.bc ~ sample.size, type="b", pch=pch_vec[2], cex=1.8, col=col_vec[2])
  lines(sim.bias.enzmann1 ~ sample.size, type="b", pch=pch_vec[3], cex=1.5, col=col_vec[3])
  lines(sim.bias.enzmann2 ~ sample.size, type="b", pch=pch_vec[4], cex=1.5, col=col_vec[4])
  lines(sim.bias.enzmann3 ~ sample.size, type="b", pch=pch_vec[5], cex=1.5, col=col_vec[5])
  lines(sim.bias.schnapp1 ~ sample.size, type="b", pch=pch_vec[6], cex=1.5, col=col_vec[6])    
  lines(sim.bias.schnapp2 ~ sample.size, type="b", pch=pch_vec[7], cex=1.5, col=col_vec[7])    
  lines(sim.bias.schnapp3 ~ sample.size, type="b", pch=pch_vec[8], cex=1.5, col=col_vec[8])    
  legend(x=1000, y=-0.01, legend=c("Naive Estimator", "Our Estimator", 
                                  "ES (No Liar)", "ES (25% Liar)", "ES (50% Liar)",
                                  "CMR-I (No Liar)", "CMR-I (25% Liar)", "CMR-I (50% Liar)"),
         col=col_vec, pch=pch_vec, box.lty=0, cex=1)


  title("Bias", line=1, cex.main=1.5)
  
  plot(sim.rmse.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, xlim=c(200,5000), ylim=c(0,0.012))
  abline(h=0, lty=2, col="gray50", lwd=2)
  lines(sim.rmse.naive ~ sample.size, type="b", cex=1.8, col=col_vec[1])
  lines(sim.rmse.bc ~ sample.size, type="b", pch=pch_vec[2], cex=1.8, col=col_vec[2])
  lines(sim.rmse.enzmann1 ~ sample.size, type="b", pch=pch_vec[3], cex=1.5, col=col_vec[3])
  lines(sim.rmse.enzmann2 ~ sample.size, type="b", pch=pch_vec[4], cex=1.5, col=col_vec[4])
  lines(sim.rmse.enzmann3 ~ sample.size, type="b", pch=pch_vec[5], cex=1.5, col=col_vec[5])
  lines(sim.rmse.schnapp1 ~ sample.size, type="b", pch=pch_vec[6], cex=1.5, col=col_vec[6])    
  lines(sim.rmse.schnapp2 ~ sample.size, type="b", pch=pch_vec[7], cex=1.5, col=col_vec[7])    
  lines(sim.rmse.schnapp3 ~ sample.size, type="b", pch=pch_vec[8], cex=1.5, col=col_vec[8])  
  title("Root Mean Square Error", line=1, cex.main=1.5)
  
  plot(sim.cover.naive ~ sample.size, type="n", ylab="", xlab="", las=1,
       main="", cex.lab=1.5, cex.main=1.8, 
       xlim=c(200,5000), ylim=c(0,100))
  abline(h=95, lty=2, col="gray50", lwd=2)
  lines(sim.cover.naive ~ sample.size, type="b", cex=1.8, col=col_vec[1])
  lines(sim.cover.bc ~ sample.size, type="b", pch=pch_vec[2], cex=1.8, col=col_vec[2])
  lines(sim.cover.enzmann1 ~ sample.size, type="b", pch=pch_vec[3], cex=1.5, col=col_vec[3])
  lines(sim.cover.enzmann2 ~ sample.size, type="b", pch=pch_vec[4], cex=1.5, col=col_vec[4])
  lines(sim.cover.enzmann3 ~ sample.size, type="b", pch=pch_vec[5], cex=1.5, col=col_vec[5])
  lines(sim.cover.schnapp1 ~ sample.size, type="b", pch=pch_vec[6], cex=1.5, col=col_vec[6])    
  lines(sim.cover.schnapp2 ~ sample.size, type="b", pch=pch_vec[7], cex=1.5, col=col_vec[7])    
  lines(sim.cover.schnapp3 ~ sample.size, type="b", pch=pch_vec[8], cex=1.5, col=col_vec[8])  
  title("Coverage of 95%\nConfidence Interval", line=1, cex.main=1.5)
  
mtext("Sample Size", side=1, line=0, outer=T, cex=1)        

dev.off()

  
######################################################
# FIGURE B2
######################################################

df <- as.data.frame(RElist.bc)
df2 <- as.data.frame(RElist.enzmann)
colnames(df) <- c("200", "500", "1000", "2000", "5000")
colnames(df2) <- c("200", "500", "1000", "2000", "5000")

pdf("FigureB2.pdf", width=11, height=4)
par(mfrow=c(1,2), mar=c(2, 2, 4.1, 1.5), oma = c(1,1,1,0) + 0.1, mgp=c(3,0.6,0)) # b l t r
boxplot(df, las=1, col="firebrick4", ylim=c(0,10))
title("Our Estimator / Naive", line=1, cex.main=1.5)
boxplot(df2, las=1, col="gray80", ylim=c(0,10))
title("Enzmann-Schnapp Correction / Naive", line=1, cex.main=1.5)
mtext(side=1, outer=T, "Sample Size")
mtext(side=3, outer=T, "Relative Length of 95% Confidence Intervals", cex=1.5, line=-1)
dev.off()


end <- Sys.time()
start - end

#############################################################################
# END OF THIS R SOURCE FILE
#############################################################################